Participant Fluctuation Amplitude Data
fluctuations.glasser.pnc <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/PNC/glasserfluctuations_demographics_finalsample.csv") #fMRI + demographics, generated with fitGAMs_fluctuationamplitude_age.RGlasser Labels
glasser.parcel.labels <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Maps/parcellations/surface/glasser360_regionlist.csv", header = T) #glasser region listSensorimotor-Association Axis
S.A.axis.cifti <- read_cifti("/cbica/projects/spatiotemp_dev_plasticity/Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo, vertex-wise axis average ranks
S.A.axis <- as.data.frame(cbind(rank(S.A.axis.cifti$data), names(S.A.axis.cifti$Parcel)))
colnames(S.A.axis) <- c("SA.axis","orig_parcelname")
S.A.axis <- merge(S.A.axis, glasser.parcel.labels, by="orig_parcelname", sort = F)
S.A.axis$SA.axis <- as.numeric(S.A.axis$SA.axis)
rm(S.A.axis.cifti)Intracortical Myelination Developmental Data
myelin.partialR2.cifti <- read_cifti("/cbica/projects/spatiotemp_dev_plasticity/Myelin/hcpd_n628_myelin_sAge_partial_bayes_r2.pscalar.nii") #T1/T2 ratio - age effect size
myelin.maxdev.cifti <- read_cifti("/cbica/projects/spatiotemp_dev_plasticity/Myelin/hcpd_n628_median_posterior_age_of_max_slope_myelination.pscalar.nii") #age of maximal T1/T2 ratio increase
myelin <- as.data.frame(cbind(myelin.partialR2.cifti$data, myelin.maxdev.cifti$data))
colnames(myelin) <- c("myelin.partialR2","myelin.maxdev")
myelin$label <- glasser.parcel.labels$label
rm(myelin.partialR2.cifti)
rm(myelin.maxdev.cifti)SNR Mask
SNR.mask <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Maps/parcellations/surface/SNRmask_glasser360.csv")Spin Test Parcel Rotation Matrix
source("/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/R/rotate.parcellation.R")
source("/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/R/perm.sphere.p.R")
glasser.coords <- read.table("/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/sphere_HCP.txt", header=F) #coordinates of glasser parcel centroids on the freesurfer sphere
perm.id.full <- rotate.parcellation(coord.l = as.matrix(glasser.coords[1:180,]), coord.r = as.matrix(glasser.coords[181:360,]), nrot = 10000) #rotate the glasser parcellation 10,000 times on the freesurfer sphere to generate spatial nulls for spin-based permutation significance testing
saveRDS(perm.id.full, "/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds")perm.id.full <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spinsGAM Results
gam.age.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_age_statistics_glasser.csv") #GAM age smooth statistics, generated with fitGAMs_fluctuationamplitude_age.R
#create df.dev with GAM statistics, S-A axis, plasticity measures, and SNR mask
df.list <- list(gam.age.glasser, S.A.axis, myelin, SNR.mask) #dfs to merge
df.dev <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list)
df.dev <- df.dev %>% select(label, orig_parcelname, everything()) #order columns
cols = c(3:15)
df.dev[,cols] = apply(df.dev[,cols], 2, function(x) as.numeric(as.character(x))) #format numerics
#create df.dev.spin formatted for spatial permutation testing
df.dev.spin <- rbind(df.dev[181:360,], df.dev[1:180,]) #format df as left hemisphere -> right hemisphere for spin tests
df.dev.spin$GAM.age.partialR2[df.dev.spin$SNR.mask == 0] <- NA #format data for spin tests by assigning NA to low SNR parcels, treating them like the medial wall
df.dev.spin$minage.decrease[df.dev.spin$SNR.mask == 0] <- NA
df.dev <- df.dev %>% filter(SNR.mask != 0) #include only high SNR parcels (N=336) in analyses
df.dev$SA.axis.bin <- as.numeric(cut2(df.dev$SA.axis, g=10)) #divide the S-A axis into 10 ranked binsgam.fitted.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_age_predictedfits_glasser.csv") #GAM predicted fits, generated with fitGAMs_fluctuationamplitude_age.R
gam.fitted.glasser <- inner_join(gam.fitted.glasser, df.dev, by="label", sort = F)gam.smoothestimates.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_age_smoothestimates_glasser.csv") #GAM smooth estimates, generated with fitGAMs_fluctuationamplitude_age.R
gam.smoothestimates.glasser <- inner_join(gam.smoothestimates.glasser, df.dev, by="label", sort = F)gam.derivatives.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_age_derivatives_glasser.csv") #GAM true model age smooth derivatives, generated with fitGAMs_fluctuationamplitude_age.R
gam.derivatives.glasser <- left_join(gam.derivatives.glasser, S.A.axis, by="label", sort = F)
gam.derivatives.glasser <- left_join(gam.derivatives.glasser, SNR.mask, by="label", sort = F)
gam.derivatives.glasser <- gam.derivatives.glasser %>% filter(SNR.mask != 0) gam.env.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_environment_statistics_glasser.csv") #GAM environment statistics, generated with fitGAMs_fluctuationamplitude_environment.R
#create df.env with GAM results, S-A axis, and SNR mask
df.list <- list(gam.env.glasser, S.A.axis, SNR.mask) #dfs to merge
df.env <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list)
#create df.env.spin formatted for spatial permutation testing
df.env.spin <- rbind(df.env[181:360,], df.env[1:180,]) #format df as left hemisphere -> right hemisphere for spin tests
df.env.spin$GAM.env.tvalue[df.env.spin$SNR.mask == 0] <- NA #format data for spin tests by assigning NA to low SNR parcels, treating them like the medial wall
df.env.spin$GAM.env.edu.tvalue[df.env.spin$SNR.mask == 0] <- NA
df.env <- df.env %>% filter(SNR.mask != 0) #include only high SNR parcels (N=336) in analysesgam.smoothestimates.glasser.lh <- gam.smoothestimates.glasser[33401:67200,]
ggplot(gam.smoothestimates.glasser.lh,aes(age,est,group=index,color=GAM.age.partialR2)) +
geom_line(size=.8, alpha = .8) +
paletteer::scale_color_paletteer_c("pals::ocean.matter", direction = -1, limits = c(-0.07, .03), oob = squish) +
theme_classic() +
labs(x = "\nAge", y = "Fluctuation Amplitude (zero centered)\n" ) +
theme(legend.position = "none") +
theme(axis.text = element_text(size=15, family = "Arial", color = c("black")), axis.title = element_text(size=15, family = "Arial", color = c("black"))) +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) Region V1
V1.smooth <- gam.fitted.glasser %>% filter(label == "rh_R_V1")
V1.derivatives <- gam.derivatives.glasser %>% filter(label == "rh_R_V1")
ggplot(data = fluctuations.glasser.pnc, aes(x = age, y = rh_R_V1)) +
geom_point(color="#FBDC9D") +
geom_line(data = V1.smooth, aes(x = age, y = fitted), color="#F8B57B", size=2) +
geom_ribbon(data = V1.smooth, aes(x = age, y = fitted, ymin = lower, ymax = upper) ,alpha = .7, linetype = 0, fill="#F8B87D",) +
labs(x="\nAge", y="Fluctuation Amplitude\n") +
theme_classic() +
theme(
axis.text = element_text(size=15, family = "Arial", color = c("black")),
axis.title.x=element_text(size=15, family ="Arial", color = c("black")),
axis.title.y=element_text(size=15, family ="Arial", color = c("black"))) +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
scale_y_continuous(breaks=c(0.5, 1, 1.5, 2, 2.5, 3))ggplot(data = V1.derivatives) +
geom_tile(aes(x = age, y = .1, fill = significant.derivative)) +
scale_fill_gradient(low = "#EFAF77", high = "#FFE0A1", limits = c(-0.0323371,-.01), na.value = "white") +
labs(x="\nAge", fill = "Derivative") +
theme_classic() +
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.line = element_blank(),
axis.ticks.y = element_blank(),
axis.title.x=element_text(size=15, family ="Arial", color = c("black")))Region p24pr
p24pr.smooth <- gam.fitted.glasser %>%filter(label == "rh_R_p24pr")
p24pr.derivatives <- gam.derivatives.glasser %>% filter(label == "rh_R_p24pr")
ggplot(data = fluctuations.glasser.pnc, aes(x = age, y = rh_R_p24pr)) +
geom_point(color="#DF5E54") +
geom_line(data = p24pr.smooth, aes(x = age, y = fitted), color="#D54D55", size=2) +
geom_ribbon(data = p24pr.smooth, aes(x = age, y = fitted, ymin = lower, ymax = upper) ,alpha = .7, linetype = 0, fill="#D54D55",) +
labs(x="\nAge", y="Fluctuation Amplitude\n") +
theme_classic() +
theme(
axis.text = element_text(size=15, family = "Arial", color = c("black")),
axis.title.x=element_text(size=15, family ="Arial", color = c("black")),
axis.title.y=element_text(size=15, family ="Arial", color = c("black"))) +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
scale_y_continuous(breaks=c(0.5, 1, 1.5, 2, 2.5, 3))ggplot(data = p24pr.derivatives) +
geom_tile(aes(x = age, y = .5, fill = significant.derivative)) +
scale_fill_gradient(high = alpha("#DF5E54",0.7), low = "#D54D55", na.value = "white", limits = c(min(p24pr.derivatives$significant.derivative),-.001)) +
labs(x="\nAge", fill = "Derivative") +
theme_classic() +
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.line = element_blank(),
axis.ticks.y = element_blank(),
axis.title.x=element_text(size=15, family ="Arial", color = c("black")))Region IFSa
#rh_R_9m
IFSa.smooth <- gam.fitted.glasser %>%filter(label == "lh_L_IFSa")
IFSa.derivatives <- gam.derivatives.glasser %>% filter(label == "lh_L_IFSa")
ggplot(data = fluctuations.glasser.pnc, aes(x = age, y = lh_L_IFSa)) +
geom_point(color="#6E195E") +
geom_line(data = IFSa.smooth, aes(x = age, y = fitted), color="#381043", size=2) +
geom_ribbon(data = IFSa.smooth, aes(x = age, y = fitted, ymin = lower, ymax = upper) ,alpha = .7, linetype = 0, fill="#381043",) +
labs(x="\nAge", y="Fluctuation Amplitude\n") +
theme_classic() +
theme(
axis.text = element_text(size=15, family = "Arial", color = c("black")),
axis.title.x=element_text(size=15, family ="Arial", color = c("black")),
axis.title.y=element_text(size=15, family ="Arial", color = c("black"))) +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
scale_y_continuous(breaks=c(0.5, 1, 1.5, 2, 2.5, 3))ggplot(data = IFSa.derivatives) +
geom_tile(aes(x = age, y = .5, fill = significant.derivative)) +
scale_fill_gradient2(high = "#6E195E", low = "#cf93c4", midpoint = 0, mid = "white", na.value = "white") +
labs(x="\nAge", fill = "Derivative") +
theme_classic() +
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.line = element_blank(),
axis.ticks.y = element_blank(),
axis.title.x=element_text(size=15, family ="Arial", color = c("black")))Fluctuation amplitude age effects map
ggseg(.data = df.dev, atlas = "glasser", mapping=aes(fill=GAM.age.partialR2), position = c("stacked"), hemisphere = c("right")) +
theme_void() +
labs(fill="Age Effect") +
paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(-0.07, .03), oob = squish)T1/T2 ratio age effects map
ggseg(.data = df.dev, atlas = "glasser", mapping=aes(fill=myelin.partialR2), position = c("stacked"), hemisphere = "right") +
theme_void() +
labs(fill="Age Effect") +
paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = 1, limits = c(0.0,0.3), oob = squish) ## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi side region label orig_parcelname GAM.age.Fvalue
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L_V1 L_V1_ROI 27.6
## 2 <NA> <NA> <NA> <NA> <NA> lh_L_MST L_MST_ROI 14.7
## 3 <NA> <NA> <NA> <NA> <NA> lh_L_V6 L_V6_ROI 7.44
## 4 <NA> <NA> <NA> <NA> <NA> lh_L_V2 L_V2_ROI 33.4
## 5 <NA> <NA> <NA> <NA> <NA> lh_L_V3 L_V3_ROI 27.8
## 6 <NA> <NA> <NA> <NA> <NA> lh_L_V4 L_V4_ROI 45.6
## 7 <NA> <NA> <NA> <NA> <NA> lh_L_V8 L_V8_ROI 29.3
## 8 <NA> <NA> <NA> <NA> <NA> lh_L_4 L_4_ROI 20.6
## 9 <NA> <NA> <NA> <NA> <NA> lh_L_3b L_3b_ROI 16.2
## 10 <NA> <NA> <NA> <NA> <NA> lh_L_FEF L_FEF_ROI 13.3
## # … with 159 more rows, and 13 more variables: GAM.age.pvalue <dbl>,
## # GAM.age.partialR2 <dbl>, Anova.age.pvalue <dbl>, age.onsetchange <dbl>,
## # age.peakchange <dbl>, minage.decrease <dbl>, maxage.increase <dbl>,
## # age.maturation <dbl>, SA.axis <dbl>, myelin.partialR2 <dbl>,
## # myelin.maxdev <dbl>, SNR.mask <dbl>, SA.axis.bin <dbl>
Correlation (rho + p.spin) between fluctuation amplitude age effects and T1/T2 ratio age effects
##
## Spearman's rank correlation rho
##
## data: df.dev$GAM.age.partialR2 and df.dev$myelin.partialR2
## S = 10557746, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.6699693
perm.sphere.p(df.dev.spin$GAM.age.partialR2, df.dev.spin$myelin.partialR2, perm.id.full,corr.type='spearman')## [1] 0.00045
Correlation plot
ggplot(df.dev, aes(x=myelin.partialR2, y=GAM.age.partialR2, fill = myelin.partialR2)) +
geom_point(color = "white",shape=21, size=3.7) +
paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = 1, limits = c(0.0,0.3), oob = squish) +
labs(x="\nMyelin Age Effect", y="Fluctuation Amplitude Age Effect\n") +
geom_smooth(method='lm', se=TRUE, fill=alpha(c("gray70"),.7), col="black") +
theme(
axis.title.x=element_text(size=15, family ="Arial", color = c("black")),
axis.title.y=element_text(size=15, family ="Arial", color = c("black")),
axis.line = element_line(color = "black"),
axis.text=element_text(size=15, family ="Arial", color = c("black")),
panel.background=element_blank(),
legend.position = "none")## `geom_smooth()` using formula 'y ~ x'
Fluctuation amplitude age of decrease onset
ggseg(.data = df.dev, atlas = "glasser", mapping=aes(fill=minage.decrease), position = c("stacked"), hemisphere = "right") +
theme_void() +
labs(fill="Age Decrease Onset") +
paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(9, 17.5), oob = squish) T1/T2 ratio age of maximal increase
ggseg(.data = df.dev, atlas = "glasser", mapping=aes(fill=myelin.maxdev), position = c("stacked"), hemisphere = "right") +
theme_void() +
labs(fill="Age Maximal Increase") +
paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(10,17.5), oob = squish) ## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi side region label orig_parcelname GAM.age.Fvalue
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L_V1 L_V1_ROI 27.6
## 2 <NA> <NA> <NA> <NA> <NA> lh_L_MST L_MST_ROI 14.7
## 3 <NA> <NA> <NA> <NA> <NA> lh_L_V6 L_V6_ROI 7.44
## 4 <NA> <NA> <NA> <NA> <NA> lh_L_V2 L_V2_ROI 33.4
## 5 <NA> <NA> <NA> <NA> <NA> lh_L_V3 L_V3_ROI 27.8
## 6 <NA> <NA> <NA> <NA> <NA> lh_L_V4 L_V4_ROI 45.6
## 7 <NA> <NA> <NA> <NA> <NA> lh_L_V8 L_V8_ROI 29.3
## 8 <NA> <NA> <NA> <NA> <NA> lh_L_4 L_4_ROI 20.6
## 9 <NA> <NA> <NA> <NA> <NA> lh_L_3b L_3b_ROI 16.2
## 10 <NA> <NA> <NA> <NA> <NA> lh_L_FEF L_FEF_ROI 13.3
## # … with 159 more rows, and 13 more variables: GAM.age.pvalue <dbl>,
## # GAM.age.partialR2 <dbl>, Anova.age.pvalue <dbl>, age.onsetchange <dbl>,
## # age.peakchange <dbl>, minage.decrease <dbl>, maxage.increase <dbl>,
## # age.maturation <dbl>, SA.axis <dbl>, myelin.partialR2 <dbl>,
## # myelin.maxdev <dbl>, SNR.mask <dbl>, SA.axis.bin <dbl>
Correlation (rho + p.spin) between fluctuation amplitude age of decrease onset and T1/T2 ratio max age of increase
## Warning in cor.test.default(df.dev$minage.decrease, df.dev$myelin.maxdev, :
## Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: df.dev$minage.decrease and df.dev$myelin.maxdev
## S = 1677077, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6382729
perm.sphere.p(df.dev.spin$minage.decrease, df.dev.spin$myelin.maxdev, perm.id.full,corr.type='spearman')## [1] 0.00125
#Age of decrease onset > 8.5
late.decrease <- df.dev %>% filter(minage.decrease > 8.5)
cor.test(late.decrease$minage.decrease, late.decrease$myelin.maxdev, method=c("spearman"))## Warning in cor.test.default(late.decrease$minage.decrease,
## late.decrease$myelin.maxdev, : Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: late.decrease$minage.decrease and late.decrease$myelin.maxdev
## S = 349670, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6402458
late.decrease.spin <- df.dev.spin
late.decrease.spin$minage.decrease[late.decrease.spin$minage.decrease < 8.5] <- NA
perm.sphere.p(late.decrease.spin$minage.decrease, late.decrease.spin$myelin.maxdev, perm.id.full,corr.type='spearman')## [1] 0.01565
Correlation plot
ggplot(df.dev, aes(x=myelin.maxdev, y=minage.decrease, fill = myelin.maxdev)) +
geom_point(color = "white",shape=21, size=3.7) +
paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(10,17.5), oob = squish) +
labs(x="\nMyelin Age of Maximal Increase", y="Fluctuation Amplitude Age of Decrease Onset\n") +
geom_smooth(method='lm', se=TRUE, fill=alpha(c("gray70"),.7), col="black") +
theme(
axis.title.x=element_text(size=15, family ="Arial", color = c("black")),
axis.title.y=element_text(size=15, family ="Arial", color = c("black")),
axis.line = element_line(color = "black"),
axis.text=element_text(size=15, family ="Arial", color = c("black")),
panel.background=element_blank(),
legend.position = "none")ggplot(late.decrease, aes(x=myelin.maxdev, y=minage.decrease, fill = myelin.maxdev)) +
geom_point(color = "white",shape=21, size=3.7) +
paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(10,17.5), oob = squish) +
labs(x="\nMyelin Age of Maximal Increase", y="Fluctuation Amplitude Age of Decrease Onset\n") +
geom_smooth(method='lm', se=TRUE, fill=alpha(c("gray70"),.7), col="black") +
theme(
axis.title.x=element_text(size=15, family ="Arial", color = c("black")),
axis.title.y=element_text(size=15, family ="Arial", color = c("black")),
axis.line = element_line(color = "black"),
axis.text=element_text(size=15, family ="Arial", color = c("black")),
panel.background=element_blank(),
legend.position = "none")Sensorimotor-association axis
ggseg(.data = df.dev, atlas = "glasser", mapping=aes(fill=SA.axis), position = c("stacked")) +
theme_void() +
theme(legend.position = "none") +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 180)Correlation (rho + p.spin) between fluctuation amplitude age effects and the sensorimotor-association axis
##
## Spearman's rank correlation rho
##
## data: df.dev$GAM.age.partialR2 and df.dev$SA.axis
## S = 2934378, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.5358554
perm.sphere.p(df.dev.spin$GAM.age.partialR2, df.dev.spin$SA.axis, perm.id.full,corr.type='spearman')## [1] 0.00215
Correlation plot
ggplot(df.dev, aes(x=SA.axis, y=GAM.age.partialR2, fill = SA.axis)) +
geom_point(color = "white",shape=21, size=4.5) +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 180) +
labs(x="\nSensorimtor-Association Axis Rank", y="Fluctuation Amplitude Age Effect\n") +
geom_smooth(method='lm', se=TRUE, fill=alpha(c("gray70"),.7), col="black") +
theme(
axis.title.x=element_text(size=15, family ="Arial", color = c("black")),
axis.title.y=element_text(size=15, family ="Arial", color = c("black")),
axis.line = element_line(color = "black"),
axis.text=element_text(size=15, family ="Arial", color = c("black")),
panel.background=element_blank(),
legend.position = "none")Correlation (rho + p.spin) between fluctuation amplitude age of decrease onset and the sensorimotor-association axis
## Warning in cor.test.default(df.dev$minage.decrease, df.dev$SA.axis, method =
## c("spearman")): Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: df.dev$minage.decrease and df.dev$SA.axis
## S = 2477911, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4655417
## [1] 0.0388
## Warning in cor.test.default(late.decrease$minage.decrease,
## late.decrease$SA.axis, : Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: late.decrease$minage.decrease and late.decrease$SA.axis
## S = 314027, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6769174
perm.sphere.p(late.decrease.spin$minage.decrease, late.decrease.spin$SA.axis, perm.id.full,corr.type='spearman')## [1] 0.0011
Correlation plot
ggplot(df.dev, aes(x=SA.axis, y=minage.decrease, fill = SA.axis)) +
geom_point(color = "white",shape=21, size=3.7) +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 180) +
labs(x="\nSensorimtor-Association Axis", y="Fluctuation Amplitude Age of Decrease Onset\n") +
geom_smooth(method='lm', se=TRUE, fill=alpha(c("gray70"),.7), col="black") +
theme(
axis.title.x=element_text(size=15, family ="Arial", color = c("black")),
axis.title.y=element_text(size=15, family ="Arial", color = c("black")),
axis.line = element_line(color = "black"),
axis.text=element_text(size=15, family ="Arial", color = c("black")),
panel.background=element_blank(),
legend.position = "none")ggplot(late.decrease, aes(x=SA.axis, y=minage.decrease, fill = SA.axis)) +
geom_point(color = "white",shape=21, size=3.7) +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 180) +
labs(x="\nSensorimtor-Association Axis", y="Fluctuation Amplitude Age of Decrease Onset\n") +
geom_smooth(method='lm', se=TRUE, fill=alpha(c("gray70"),.7), col="black") +
theme(
axis.title.x=element_text(size=15, family ="Arial", color = c("black")),
axis.title.y=element_text(size=15, family ="Arial", color = c("black")),
axis.line = element_line(color = "black"),
axis.text=element_text(size=15, family ="Arial", color = c("black")),
panel.background=element_blank(),
legend.position = "none")Sensorimotor-association axis binned GAM smooths
gam.smoothestimates.glasser$age <- round(gam.smoothestimates.glasser$age, 3)
meansmooth.bins <- gam.smoothestimates.glasser %>% group_by(age, SA.axis.bin) %>% do(est.mean = mean(.$est)) %>% unnest(cols = c(est.mean)) #calculate the average smooth estimate at each age within each of the 10 S-A axis binsggplot(data = meansmooth.bins, aes(x = age, y = est.mean, group = SA.axis.bin, color = as.factor(SA.axis.bin))) +
geom_line(size = 2.5) +
scale_color_manual(values = c("#FFC228", "#FFCA4E", "#FFD16A", "#FFDA89", "#FFE6B0", "#D7BCDA", "#BE93C4", "#9859A4", "#863E95", "#6F1382")) +
theme_classic() +
labs(x = "\nAge", y = "Fluctuation Amplitude (zero centered)\n", color = "Sensorimotor-Association Axis Bin") +
theme(axis.text = element_text(size=15, family = "Arial", color = c("black")), axis.title = element_text(size=15, family = "Arial", color = c("black"))) +
theme(legend.position = c(.585, 0.94), legend.direction = "horizontal", legend.text = element_text(size=13), legend.title = element_text(size=12, family = "Arial", color = c("black"))) +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))Principal component analysis to identify the spatial axis that captures the most variance in regional fluctuation amplitude developmental trajectories
gam.smoothestimates.glasser <- gam.smoothestimates.glasser %>% select(label, age, est)
gam.smoothestimates.long <- gam.smoothestimates.glasser %>% pivot_wider(names_from="label",values_from="est")
gam.smoothestimates.long <- gam.smoothestimates.long %>% select(-age) #ages by regions matrix for PCA
smooths.pca <- prcomp(gam.smoothestimates.long) #PCAVariance explained by PC1
## [1] 0.87026
Developmental PC1 correlation with the sensorimotor-association axis
PC1 <- as.data.frame(smooths.pca$rotation[,1]) #first principal component loadings
PC1$label <- row.names(PC1)
colnames(PC1) <- c("PC1","label")
PC1$ranked <- rank(PC1$PC1)
cor.test(PC1$PC1, df.dev$SA.axis, method=c("spearman"))##
## Spearman's rank correlation rho
##
## data: PC1$PC1 and df.dev$SA.axis
## S = 1896048, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.700093
df.dev.spin <- left_join(df.dev.spin, PC1, by="label")
perm.sphere.p(df.dev.spin$PC1, df.dev.spin$SA.axis, perm.id.full,corr.type='spearman')## [1] 0
Correlation plot
PC1 <- left_join(PC1, S.A.axis, by="label")
ggplot(PC1, aes(x=SA.axis, y=PC1, fill = SA.axis)) +
geom_point(color = "white",shape=21, size=3.7) +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 180) +
labs(x="\nSensorimtor-Association Axis", y="Development PC1\n") +
geom_smooth(method='lm', se=TRUE, fill=alpha(c("gray70"),.7), col="black") +
theme(
axis.title.x=element_text(size=15, family ="Arial", color = c("black")),
axis.title.y=element_text(size=15, family ="Arial", color = c("black")),
axis.line = element_line(color = "black"),
axis.text=element_text(size=15, family ="Arial", color = c("black")),
panel.background=element_blank(),
legend.position = "none")PCA brain map
ggseg(.data = PC1, atlas = "glasser", mapping=aes(fill=PC1), position = c("stacked")) + theme_void() +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.055, limits = c(-0.08,0.008), oob = squish) + theme(legend.position = "none")## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi side region label PC1 ranked orig_parcelname SA.axis
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <chr> <dbl>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L_10pp 0.0175 319 L_10pp_ROI 282
Regional derivative by age plot
gam.derivatives.glasser.lh <- gam.derivatives.glasser[33401:67200,]
ggplot(gam.derivatives.glasser.lh,aes(age,derivative,group=index,color=SA.axis)) +
geom_line(size=1, alpha = .7) +
scale_color_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, midpoint = 180) +
theme_classic() +
theme(legend.position = "none") +
labs(x = "\nAge", y = "Derivative\n", color = "Sensorimotor-Association Axis Bin") +
theme(axis.text = element_text(size=15, family = "Arial", color = c("black")), axis.title = element_text(size=15, family = "Arial", color = c("black"))) +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))Age 10 derivative map
age10.derivs <- gam.derivatives.glasser %>% filter(age == 10.03015)
age10.derivs$derivative <- rank(age10.derivs$derivative)
ggseg(.data = age10.derivs, atlas = "glasser", mapping=aes(fill=derivative), position = c("stacked"), hemisphere = "right") +
theme_void() +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 168)Age 15 derivative map
age15.derivs <- gam.derivatives.glasser %>% filter(age == 15.02429)
age15.derivs$derivative <- rank(age15.derivs$derivative)
ggseg(.data = age15.derivs, atlas = "glasser", mapping=aes(fill=derivative), position = c("stacked"), hemisphere = "right") +
theme_void() +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 168)Age 20 derivative map
age20.derivs <- gam.derivatives.glasser %>% filter(age == 20.01843)
age20.derivs$derivative <- rank(age20.derivs$derivative)
ggseg(.data = age20.derivs, atlas = "glasser", mapping=aes(fill=derivative), position = c("stacked"), hemisphere = "right") +
theme_void() +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 168)Derivative by sensorimotor-association axis plot
ggplot(data=gam.derivatives.glasser,aes(x = age,y = SA.axis, group = index)) +
geom_line(aes(color=derivative), size=1) +
#paletteer::scale_color_paletteer_c("grDevices::RdYlBu", direction = -1, limits=c(-0.03,0.03), oob = squish) +
paletteer::scale_color_paletteer_c("grDevices::Sunset", direction = -1, limits=c(-0.03,0.03), oob = squish) +
#paletteer::scale_color_paletteer_c("grDevices::Inferno", direction = -1, limits=c(-0.03,0.03), oob = squish) +
theme_classic() +
ylab("Sensorimotor-Association Axis Rank\n") +
xlab("\nAge") +
labs(color="Amplitude change per year\n") +
theme(axis.text.x = element_text(size = 15, color = c("black"))) +
theme(axis.text.y = element_text(size = 15, color = c("black"))) +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))Posterior derivative - sensorimotor-association axis correlations, N=10,000
deriv.SAaxis.posteriorcorrs <- read.table("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/SAaxis_posteriorderivative_correlation_byage_glasser.csv", header = F, sep = ",") #get the regional derivative-regional axis rank correlation at each age for all draws from the posterior distribution
colnames(deriv.SAaxis.posteriorcorrs) <- sprintf("draw%s",seq(from = 1, to = ncol(deriv.SAaxis.posteriorcorrs)))
deriv.SAaxis.posteriorcorrs <- cbind(deriv.SAaxis.posteriorcorrs, (gam.derivatives.glasser %>% group_by(age) %>% group_keys()))age.max.corr <- deriv.SAaxis.posteriorcorrs %>% #find the age at which the correlation is largest for each draw
summarise(across(contains("draw"),
.fns = function(x){
round(deriv.SAaxis.posteriorcorrs$age[which.max(x)],4)
}))
age.max.corr <- t(age.max.corr)
age.max.corr.median <- median(age.max.corr) #median age #bayes
age.max.corr.CI <- quantile(age.max.corr, probs = c(0.025, 0.975)) #compute the credible interval based on all draws
age.max.corr.lower <- age.max.corr.CI[[1]]
age.max.corr.upper <- age.max.corr.CI[[2]]## [1] 15.0243
## 2.5% 97.5%
## 14.6516 15.3224
deriv.SAaxis.posteriorcorrs <- deriv.SAaxis.posteriorcorrs %>% select(contains("draw"))
deriv.SAaxis.mediancorr <- apply(deriv.SAaxis.posteriorcorrs, 1, function(x){median(x)}) #median correlation value at each age
cor.credible.intervals <- apply(deriv.SAaxis.posteriorcorrs, 1, function(x){quantile(x, probs = c(0.025, 0.975))}) #compute the credible interval for the correlation value at each age based on all draws
cor.credible.intervals <- t(cor.credible.intervals)
cor.credible.intervals <- as.data.frame(cor.credible.intervals)
cor.credible.intervals <- cbind(cor.credible.intervals, (gam.derivatives.glasser %>% group_by(age) %>% group_keys()))
cor.credible.intervals$SA.correlation <- deriv.SAaxis.mediancorr
colnames(cor.credible.intervals) <- c("lower","upper","age","SA.correlation")cor.credible.intervals$max.corr.CI <- (cor.credible.intervals$age > age.max.corr.lower & cor.credible.intervals$age < age.max.corr.upper) #add a column that indicates whether each age is included in the age of maximal correlation credible interval (T/F)
cor.credible.intervals$max.cor.window <- cor.credible.intervals$age*cor.credible.intervals$max.corr.CI #add a column that only includes ages in this interval
cor.credible.intervals$max.cor.window[cor.credible.intervals$max.cor.window == 0] <- NA
cor.credible.intervals$zero.corr.CI <- (cor.credible.intervals$lower < 0 & cor.credible.intervals$upper > 0) #add a column that indicates whether the credible interval for the correlation includes 0
cor.credible.intervals$zero.corr.window <- cor.credible.intervals$age*cor.credible.intervals$zero.corr.CI #add a column that only includes ages in the zero window
cor.credible.intervals$zero.corr.window[cor.credible.intervals$zero.corr.window == 0] <- NA## [1] 19.12395
## [1] 22.47822
Sensorimotor-association axis development sliding window plot
ggplot(cor.credible.intervals, aes(x = age, y = SA.correlation, ymin = lower, ymax = upper)) +
geom_line(size=1.5) +
geom_ribbon(alpha=.2, fill=c("grey30")) +
labs(x="\nAge", y="Developmental Alignment the Axis\n") +
geom_ribbon(aes(x = max.cor.window, y = SA.correlation), fill="#F4A674") +
geom_ribbon(aes(x = zero.corr.window, y = SA.correlation), fill="#952162") +
theme(
axis.title.x=element_text(size=15, color = "black"),
axis.title.y=element_text(size=15, color = "black"),
axis.line = element_line(color = "black"),
axis.text=element_text(size=15, color = "black"),
panel.background=element_blank(),legend.position = "none") +
scale_y_continuous(breaks=c(0,0.2,0.4,0.6)) +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), limits = c(8, 23), expand = c(0,0))Number of cortical regions showing significant environment effects
## [1] 141
## [1] 106
Percent of cortical regions showing significant environment effects
## [1] 41.96429
## [1] 31.54762
Across-cortex environment effect size and direction
ggseg(.data = df.env, atlas = "glasser", mapping=aes(fill=GAM.env.tvalue), position = c("stacked")) +
theme_void() +
labs(fill="Age Effect") +
scale_fill_gradient2(low= "#F4A674", mid = "white", high = "#952162", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0, limits=c(-4,4), oob=squish)##
## Spearman's rank correlation rho
##
## data: df.env$GAM.env.tvalue and df.env$SA.axis
## S = 3274572, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4820453
## [1] 0
#significant regions
df.env.spin$GAM.env.tvalue.significant <- df.env.spin$GAM.env.tvalue*(p.adjust(df.env.spin$Anova.env.pvalue, method=c("fdr")) < 0.05)
df.env.spin$GAM.env.tvalue.significant[df.env.spin$GAM.env.tvalue.significant == 0] <- NA
cor.test(df.env.spin$GAM.env.tvalue.significant, df.env.spin$SA.axis, method=c("spearman"))##
## Spearman's rank correlation rho
##
## data: df.env.spin$GAM.env.tvalue.significant and df.env.spin$SA.axis
## S = 204474, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.5623229
perm.sphere.p(df.env.spin$GAM.env.tvalue.significant, df.env.spin$SA.axis, perm.id.full, corr.type='spearman')## [1] 0
Correlation plot
ggplot(df.env, aes(x=SA.axis, y=GAM.env.tvalue, fill = SA.axis)) +
geom_point(color = "white",shape=21, size=3.7) +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 180) +
labs(x="\nSensorimtor-Association Axis", y="Environment Effect\n") +
geom_smooth(method='lm', se=TRUE, fill=alpha(c("gray70"),.7), col="black") +
theme(
axis.title.x=element_text(size=15, family ="Arial", color = c("black")),
axis.title.y=element_text(size=15, family ="Arial", color = c("black")),
axis.line = element_line(color = "black"),
axis.text=element_text(size=15, family ="Arial", color = c("black")),
panel.background=element_blank(),
legend.position = "none")##
## Spearman's rank correlation rho
##
## data: df.env$GAM.env.edu.tvalue and df.env$SA.axis
## S = 2914840, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.5389458
perm.sphere.p(df.env.spin$GAM.env.edu.tvalue, df.env.spin$SA.axis, perm.id.full, corr.type='spearman')## [1] 0
#significant regions
df.env.spin$GAM.env.edu.tvalue.significant <- df.env.spin$GAM.env.edu.tvalue*(p.adjust(df.env.spin$Anova.env.edu.pvalue, method=c("fdr")) < 0.05)
df.env.spin$GAM.env.edu.tvalue.significant[df.env.spin$GAM.env.edu.tvalue.significant == 0] <- NA
cor.test(df.env.spin$GAM.env.edu.tvalue.significant, df.env.spin$SA.axis, method=c("spearman"))##
## Spearman's rank correlation rho
##
## data: df.env.spin$GAM.env.edu.tvalue.significant and df.env.spin$SA.axis
## S = 83236, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.5806434
perm.sphere.p(df.env.spin$GAM.env.edu.tvalue.significant, df.env.spin$SA.axis, perm.id.full, corr.type='spearman')## [1] 0
Correlation plot
ggplot(df.env, aes(x=SA.axis, y=GAM.env.edu.tvalue, fill = SA.axis)) +
geom_point(color = "white",shape=21, size=3.7) +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 180) +
labs(x="\nSensorimtor-Association Axis", y="Environment Effect\n") +
geom_smooth(method='lm', se=TRUE, fill=alpha(c("gray70"),.7), col="black") +
theme(
axis.title.x=element_text(size=15, family ="Arial", color = c("black")),
axis.title.y=element_text(size=15, family ="Arial", color = c("black")),
axis.line = element_line(color = "black"),
axis.text=element_text(size=15, family ="Arial", color = c("black")),
panel.background=element_blank(),
legend.position = "none")